Lya Transfer in a thick, dusty, and static medium 

1 Sang-Hyeon Ahn, 2 Hee-Won Lee, 1 Hyung-Mok Lee 
1 Astronomy Program, SEES, Seoul National University, Seoul, Korea 
2 Department of Astronomy, Yonsei University 



Received 13 June 2000; accepted 



Submitted to JKAS 2000 



- 2 - 
ABSTRACT 

We developed a Monte Carlo code that describes the resonant Lya line transfer in 
an optically thick, dusty, and static medium. The code was tested against the analytic 
formula derived by Neufeld (1990). We explain the line transfer mechanism for a wide 
range of line center optical depths by tracing histories of photons in the medium. We 
find that photons escape from the medium by a series of wing scatterings, during 
which polarization may develop. We applied our code to examine the amount of dust 
extinction around the Lya in primeval galaxies. Brief discussions on the astrophysical 
application of our work are presented. 

Subject headings: cosmology: miscellaneous — galaxies: formation — polarization - 
radiative transfer 
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1. Introduction 

Recent advances of large telescopes and CCD techniques have been unveiling very remote 
astrophysical objects such as primeval star forming galaxies and damped Lya galaxies. The 
Lya emission profiles of these primeval galaxies are categorized into three types : (1) pure Lya 
emission, (2) asymmetric or P Cyg type one, (3) broad absorption in damping wings. Moreover, 
local starburst galaxies also show similar characteristics (e.g. Kunth et al. 1998). We will not 
list up the full inventory of the starburst galaxies mentioned above, but the issue was discussed 
in our previous paper (Lee & Ahn 1999). It is plausible that dusts play an important role in the 
formation of the Lya line profile. Since Lya sources in a starburst galaxies are usually embeded in 
optically thick and dusty media, we must investigate the effects of dusts on the radiative transfer 
and line formation of Lya in the media. 

Recently Meurer et al. (1999) studied the dust extinction in local starburst galaxies, 
and applied the result to the primeval galaxies in the Hubble Deep Field. They found a tight 
relationship between the f3 indexes of ultraviolet (UV) continua and the ratios of far-infrared 
(FIR) to UV fluxes of the galaxies. The relationship is established because dusts absorb Lya line 
photons and UV continua of Ao > 912A, and then re-emit in FIR. Leitherer et al. (1999) suggested 
that the Lya luminosity is as large as 10% of the total UV luminosity in the star forming regions. 
Hence it is necessary to know what fraction of Lya photons is subject to dust extinction during 
the transfer in a very thick neutral medium. 

Studies of the Lya line transfer in an optically thick and static medium has a long history. 
Unno (1955) formulated the Lya line transfer in a dust-free medium, and Osterbrock (1962) 
proposed a simple physical picture for understanding the resonance line transfer in a thick 
medium. Adams (1972) revised Osterbrock's picture by a Monte Carlo method and gave a 
heuristic explanation on the problem. In an analytical way, Harrington (1973) solved the problem, 
for which Neufeld (1990) also gave a more general solution. However, the studies thus far are 
limited to the cases where diffusion approximation is valid and little attention has been paid to 
the polarization of the Lya flux emergent from those media with anisotropic geometry and/or 
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kinematics. More realistic treatment of the Lya line transfer is clearly necessary. 

We have developed a sophisticated Monte Carlo code that describes the radiative trasfer of 
Lya resonnance photons in optically thick, dusty, and static media, ranging from moderate optical 
depths to extremely large optical depths. In Section II, we describe the basic theory including the 
numerical method of our code. The basic results are presented in Section III, and the possible 
applications of our code are described in Section IV. 

2. Basic Theory 

2.1. Configuration and Review 

The medium used in our investigation is plane-parallel, in which hydrogen atoms and dust 
particles are mixed homogeneously. In Fig. 1, we show the configuration. We consider the case 
that a source is located in the mid-plane and radiates photons isotropically. However, it is easy to 
extend to other configurations. 

Osterbrock (1962) proposed a simple physical picture based on a simplified treatment of 
random walk processes, and predicted that (N) oc Tq , where (N) is the mean number of scattering 
and To is the line center optical depth. But Adams (1972) solved the problem by a Monte Carlo 
method and showed that (N) oc To rather than (TV) oc Tq, where a is the Voigt parameter. 
Also (N) and the peak of emergent spectra is dependent upon citq rather than tq alone. He 
made an explanation that the line transfer process is dominated by wing scattering, and is well 
approximated by random walk processes. He also pointed out that the transfer of Lya photons in 
a thick medium can be approximated as a diffusion process both in real space and in frequency 
space. He suggested that in moderately thick media (i.e. 10 3 < tq < 10 3 /a) Lya photons escape 
by a single longest flight, and for extremely thick media (i.e. tq > 10 3 /a) by a single longest 
excursion. 

Harrington (1973) solved the radiative transfer equation presented by Unno (1955) in an 
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Fig. 1. — Schematic diagram describing the configuration in our investigation. We adopt a plane- 
parallel and uniform medium of optical depth To to both sides in z-direction. The Lya photons are 
generated isotropically from the source located at the center of the medium. 
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analytic manner. He adopted the diffusion approximation and used the Sturm-Liouville theory to 
obtain the analytic solution. However, he considered only for one specific configuration, in which a 
monochromatic (Ao = 1216A) source is located at the center of a plane-parallel medium. Neufeld 
(1990) gave analytic solutions for a more generalized version of the problem, where the source 
produces arbitrary initial frequencies and various dust extinctions were included. 

In the present paper, we will consider the same problems with more general optical depths 
using Monte Carlo techniques, and the results obtained by Neufeld (1990) will be used to assure 
the validity of our method. However, it is noted that our Monte Carlo approach is more faithful 
to the problem by resorting to less approximation than other approaches mentioned above. 



2.2. Optical Depth 

In the configuration stated in the previous subsection, the optical depth is given by 



oc 



t v (s) = dl dv z n(v z )a u . (1) 



o 

In a static medium, the Lya scattering cross section is given by 

a - f T/An2 (2\ 

m e c {v — vqY + {Y / Att) z 

where m e is the electron mass, c the light velocity, /12 = 0.4126 the oscillator strength for hydrogen 

Lya, T = A 2 \ = 3.1 x 10 9 s _1 the damping constant, and vq the frequency of the line center. When 

we include the thermal motions of scatterers, v is tranformed by v' = v — fo^f, and the thermal 

motion is described by the Maxwellian distribution, 

n(v z )dv z = n HI —= — e v th dv z . (3) 

V^Vth 

Here v\ h = 2kBT g /rriH, where ran is the hydrogen mass, ks the Boltzmann constant, nn the 
neutral hydrogen density, and T g the gas temperature. 



Substituting Eqs.(2) and (3) into Eq.(l), we obtain the expression for the optical depth, 

tt^V , 1 a f°° , e- u2 

T x {s) = JuriHis / du 

m e c v D 7r 7_c 



(x — u) 2 + a 2 
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where T gA = T g /10 4 K, 



lAlxW- 13 T g4 1/2 n HI sH(x,a), (4) 



vd = —^o, (5) 
X - Au D Av D ' W 



and 



o = ■ (8) 

Here a) is the Hjerting function or the H-function, and the specific values in Rybicki & 
Lightman (1979) are used. Also the Voigt parameter a is given by 

a = 4.7 x 1(T 4 T 9 ~ 1/2 . (9) 

Therefore, the total optical depth for a given system is 

t x = 1.41 x 10- 13 r- 4 1/2 iVH/F(x, a) = T H(x,a), (10) 

where the column density of neutral hydrogen, Nhi = nniL with L being the physical size of the 
slab, and tq = 1.41 x 10 _13 T Nhi is the line center optical depth. 



2.3. Monte Carlo Code 

In this subsection, we present a detailed description of the Monte Carlo procedure. There 
are a few Monte Carlo approaches to the resonance line transfer in an optically thick and static 
medium in the literature (e.g. Adams 1972, Harrington 1973, Gould & Weinberg 1996). 

The Monte Carlo code begins with the choice of the frequency and the propagation direction 
ki of an incident photon from an assumed Lya profile, which is assumed to be monochromatic in 
our study. 



-8- 



Then we determine the next scattering site separated from the initial point by the propagation 
length defined by 



where the optical depth is assumed to be composed of hydrogen {tqH(x,o)) and dust (r<j) parts. 
Here, Td is the dust optical depth which is treated as a constant for the emission line and whose 
Galactic value is given by Draine & Lee (1984), and r = — ln(R), where R is a uniform random 
number in the interval (0, 1) generated by a subroutine ran2() suggested by Press et al.(1989). 

The emitted photon traverses a distance S, and is scattered off by hydrogen atoms until 
\z\ > ro for a slab geometry. In this scattering event the frequencies of the absorbed photon and 
the re-emitted one in the rest frame of the scatterer should be matched. 

In contrast to the case of a thin medium, for a very thick medium the scattering in the 
damping wings is not negligible. A careful treatment needs to be exercised to distinguish the 
scattering in the damping wings from the resonance scattering, because they show quite different 
behaviors in the properties including the scattering phase function and the polarization (Lee & 
Blandford 1997). 

Because the natural line width is much smaller than the Doppler width, the local velocity 
of the scatterers that can resonantly scatter the incident photon is practically a single value. 
However, when the scattering occurs in the damping wings, the local velocity of the scatterer may 
run a rather large range. Therefore, in order to enhance the efficiency of the Monte Carlo method, 
it is desirable to determine the scattering type before we determine the local velocity u of the 
scatterer. 

We present a more quantitative argument about the preceding remarks. Under the condition 
that a given photon is scattered by an atom located at a position s, the local velocity component 
u along the direction ki is chosen from the normalized distribution, 



S = r s = 



r 



T 



(11) 



H{x,a) +T d /T H(x,a) 



— w 



.2 
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e 



—H(x, a) 



(12) 



(x — u) 2 + a 2 [a 
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Here, the Hjerting function or the Voigt function H(x, a) is evaluated by a series expansion in a, 
i.e., 

H(x,a) = Ho(x,a) + aH\(x,a) + a 2 H2(x,a) 

+ a 3 H 3 (x,a) + •••, (13) 

where H n (x, a), n = 0, 1, 2, 3 are tabulated by Gray (1992). 

Because of the smallness of a, the function / has a sharp peak around u « x, for which 
the scattering is resonant. Therefore, the probability P r that a given scattering is resonant is 
approximately given by 

/oo p _a;2 
d(Au) 2 
-oo (Aw) 2 + a 2 

i?(x, a) 

The probability that scattering occurs in the damping wings is 

P nr = 1 - P r . (15) 

In the code we determine the scattering type in accordance with the scattering type 
probabilities P r and P nr . If scattering is chosen to be resonant, then we set u = x. Otherwise, the 
scattering occurs in the damping wings, and u is chosen in accordance with the velocity probability 
distribution given by Eq.(12). 

We give the propagation direction fcj of a scattered photon in accordance with the phase 
function completely faithful to the atomic physics (Lee & Ahn 1998). The scattered velocity 
component v± perpendicular to the initial direction ki on the plane spanned by ki and kf is also 
governed by the Maxwell-Boltzmann velocity distribution, which is numerically obtained using the 
subroutine gasdevQ suggested by Press et al. (1989). The contribution Ax of the perpendicular 
velocity component v± to the frequency shift along the direction of kf is obviously 

Ax = v ± [l-(k i -k f ) 2 ] 1 / 2 /c. (16) 



—H(x, a) 
a 
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Therefore, the frequency shift Xf of the scattered photon is given by 

Xf = Xi - u + u(ki ■ kf) +v±[l- {ki ■ kf) 2 } 1/2 , (17) 

where Xi is the frequency shift of the incident photon. 

In each scattering event the position of scattered photon is checked, and its path length is 
added. If the photon escapes from the medium or \z\ < To, we collect that photon according to its 
frequency and escaping direction. In collecting photons, we add the weighted fraction considering 
the dust extinction corresponding to the total path length of the photon. The whole procedure is 
repeated until typically about 10 3 photons in each frequency bin are collected. 

3. Results 

3.1. Tracing the Scattering Processes 

Adams (1972) gave a physical picture that describes the Lya line transfer in extremely thick 
media of r > 10 3 /a, where a = 4.71 x 1(T 4 for T g = 10 4 K. 

Initially a core photon becomes a wing photon after a sufficient number of scatterings and 
encountering with a violently moving atom. Due to the small optical depths, wing photons traverse 
much further in physical space than core photons. When To < 10 3 , a photon can directly escape 
from a medium once it becomes a wing photon. This process is called 'a single longest flight' by 
Adams (1972). 

As the line center optical depth gets larger to be an intermediate optical depth, 
10 3 < tq < 10 3 /a, wing photons cannot directly escape from the medium, but experience a 
large number of core scatterings followed by a few wing scatterings. However, after these wing 
scatterings, photons may return back to the core because of the so called 'the restoring force' 
(Osterbrock 1962). 

For the case of the extremely thick medium we can describe this process using a diffusion 
approximation, where the rms frequency shift per scattering is the thermal Doppler width (v t h) 
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Fig. 2. — History of a photon, x-axis is the scattering number, and we show frequencies of the 
photon in Panel A, its locations in panel B, and its type of scattering in panel C. In panel C, '0' 
means wing scattering, and '1' means core scattering. 
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Fig. 3. — History of a photon just before escape. Panels represent the same quantities to those in 
Fig. 2. We can see that the photon experiences wing scattering just before escape. 
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or x = 1, and the mean frequency shift per scattering is Ax = — l/\x\ that plays the role of the 
restoring force. This restoring force is caused by relatively large probability that photons scatter 
in the core. However, for a medium with an intermediate optical depth, the first variance of 
frequency is not quite large because the optical depth is not enough to enlarge the probability 
for wing scattering at frequencies far from the line center. Hence, photons often experience a 
small number of single longest flights before escape, and we may call these processes as 'random 
wandering'. Therefore, in this particular case, we think that a simple diffusion approximation is 
not adequate and that a more accurate approach such as a Monte Carlo method is needed. 

For a medium of an even higher optical depth, tq > 10 3 /a, a series of wing scatterings occur 
when wing photons are generated. In this case the wandering occurs both in real space and in 
frequency space. This process is called 'an excursion,' and photons escape the medium by 'a single 
longest excursion' (Adams 1972). 

These processes can be seen in our Monte Carlo results. In Fig. 2 we have shown the 
trajectory of a photon that transfers in a medium with tq = 10 4 . At the 8055th scattering, the 
photon is scattered off by an atom that moves fast and became a wing photon. At the same time, 
the photon traverses a longer distance in physical space, and also wing scattering happens. 

Before escaping from the medium, the photon experiences a series of wing scatterings. In 
Fig. 3 we show the history of the same photon, and here we can see that the last scattering 
occurs in the wing. We checked that the larger To, the larger the number of last wing scatterings, 
which may be regarded as 'a single longest excursion.' These last wing scatterings possibly induce 
polarization because wing scattering has a Rayleigh phase function (Stenflo 1980). We have 
conducted an investigation on this topic, and the results will be published in another paper. 

3.2. Emergent Lyct Profile for Dust Free Media 

Neufeld (1990) derived an analytic solution, for (aro) 1 / 3 ~> 1. According to Unno (1955), 
Harrington (1973), and Neufeld (1990), for the cases of dust-free media, the Lya line transfer 
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Fig. 4. — Our emergent profiles (solid lines) are compared with Neufeld's analytic solution (dotted 
lines), x axis is a frequency in units of the thermal width, and the total flux of the line is normalized 
to 1/87T in accordance with Neufeld's normalization. The profiles are symmertic about the origin, 
Av/Av D = 0. 



-15- 



equation is given by 



dr 2 da 2 An 

where 



d 2 J d 2 J ,E(t,x) , , 

+ ^ = -3^0-^, (18) 



'-[*) W (19) 



and <^(s) is the normalized Voigt function given by 



Here the boundary conditions are 



<p(x) = -^=H(x,a). (20) 



J(r, ±oo) = (21) 



2 /dJ\ 



J( ± r„,,) = ± 2 ff ( ± r,: C ) = ^^j ± ^ (22) 

Since the source function E(t,x) in the right hand side of Eq. (19) is very sharply peaked, 
we can apply the Green function method to the problem. With an approximaton 

(j)(x) ~ a/irx 2 (23) 

and an application of the Sturm-Lioville theory, Neufeld (1990) derived a solution for the case of a 
mid-plane source, 

J(±T °' X) = '24^cosh[(7r 4 /54) 1 /2(| x 3_ :r 3| /aro) ]- ( 24 ) 

We perform Monte Carlo calculations for a monochromatic source with Xi = 0, and compare 
the result with Eq.(24). Fig. 4 shows our results for r = 10 5 and tq = 10 6 . The solid lines are 
the results of Monte Carlo calculation, and the dotted lines are those obtained by Neufeld. The 
shapes of the profiles agree with each other, and we emphasize that our Monte Carlo code is 
written to incorporate all the quantum mechanics associated with both resonant and non-resonant 
Lya scattering. 
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The figure shows an overall disagreement: the Monte Carlo solution appears to be translated 
from the Neufeld's solution by an amount of x ~ 3, which corresponds to the frequency shift where 
the wing scattering becomes important and the profile function (p(x) can be nicely approximated 
by Eq. (24). In fact, in the cases of dust free media, no loss of Lya line photons is permitted, 
which guarantees the flux conservation. The profile function <p{x) can not be approximated by 
Eq.(23) at x < 3, and therefore, Neufeld's calculation underestimates the amount of core photons 
that are removed and ultimately redistributed to the wing regimes, which reduces the amount of 
diffusively trasferred wing photons. 



3.3. Survival Fraction for Dusty Media 



Neufeld (1990) also provided an analytic solution for the problem of the radiative transfer in 
a dusty medium. Here we will concentrate on the continuum absorption case. Harrington (1973) 
and Neufeld (1990) formulated this problem by 



d 2 J d 2 J 
Or 2 + da 2 



= -3[<f>(x) + fax)] 



E(t,x) 
Air 



P(x)J 



where 



/ 2 y/2 , 
G= \V Jo 



dx' 



(25) 



(26) 



(t>(x') + f3t{xy 

We assume that the total dust opacity (3t(x) « becasue (j)(x) » (3f Here, j3 is an absorptive part 
of dust opacity per atomic hydrogen, which is given by 

/3=1.0xlO- 8 T 9 1 4 /2 ^-^ (27) 

Here 



nHl 



(I -A) 



Nhi 



G 



(28) 



where A is albedo and the subscript G represents the Galactic value. A dust opacity per a 
hydrogen atom is given by 

Td 



Nhi 



fd(Vsca + dabs), 



(29) 
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Fig. 5. — Survival fraction of photons in our Monte Carlo calculations for dusty media of To = 10 is 
compared with previous works. Long dashed line stands for the case of a variable effective frequency, 
and dotted line for the case of a constant effective frequency. Short dashed line represents Neufeld's 
derivation for Yq by considering average scattering number of photons in media with intermediate 
optical depth. 
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where is dust optical depth, is dust-to-gas ratio, a sca and a a b s are the scattering and the 
absorption cross sections, respectively. According to Draine and Lee (1984), the mean value of the 
dust optical depth for Our Galaxy is given by Td/Nni = 1.6 x 10 -21 at Ao = 1216A 

In our Monte Carlo code, we give r^, A, as well as Nhi or tq. Therefore we can compare 
these two kinds of calculations. Neufeld also gave solutions for two kinds of approximation: (1) 
variable effective scattering frequency, (2) constant effective scattering frequency. In the former 
approximation, the survival fraction of photons is given by 



fe = 



i 



1 + SV fc! 3-7-ll----(4fc-l) 



(30) 



where y = §(f ^i/e^^/a. 

In the latter approximation, the total fraction of photons which escape the slab is given by 

1 



fe 



cosh(Yo) ' 



(31) 



where 



Yn 



3(aT y^PT 
vrC 2 



I 1/2 



(32) 



and £ is determined to be C = 0.525 by fitting the results of Hummer and Kunasz (1980) to the 
following relation between peak frequency (x s ) and the aro, 



x s = C(aro) 1 / 3 . 



(33) 



The other way of determining Yq is provided by Neufeld (see Eq. 4.35 in his paper), 
who considered the mean scattering length in the result of Hummer and Kunasz (1980). The 
consideration leads to an survival fraction 1/ cosh(Yo) with 

Y = (31.8/3r ) 1/2 . (34) 

He argued that this relationship may be applicable even in the range of intermediate optical depth, 
10 3 < r < 10 3 /o. 
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In Fig. 5 we show the main result on the survival fraction of Lya photons in a dusty medium 
with ro = 10 6 and T = 10 4 K. Our results are in good agreement with Neufeld's approximate 
solutions and those of Hummer and Kunasz. Our results are also consistent with Neufeld's 
approximation solution for the approximation of a constant effective frequency. 

On the other hand, although Neufeld noted that Eq.(34) above is valid even in intermediate 
optical depths, according to our results, this statement is not true. 

4. Discussion 

In this paper we described the Lya line transfer in an optically thick, dusty, and static 
medium. We took accurate atomic physical considerations in the Monte Carlo code, and examined 
the effect of dusts on Lya line formation in optically thick media. 

For dust free media we confirmed the line transfer mechanism, in which a number of longest 
flights occur. We emphasize that photons experience a series of wing scatterings at the moment 
of escape, at which stage significant polarization can be developed. We also confirmed that line 
profiles from our Monte Carlo calculations agree well with the analytic solution derived by Neufeld 
(1990). 

For dusty media, we calculated the survival fraction of photons for To = 10 6 and T = W 4 K, 
and found that our results are in good agreement with those obtained by Hummer and Kunasz 
(1980) and with those of Neufeld's approximation solution. 

The computational speed is limited for extremely thick media with aro > 10 3 . This difficulty 
can be overcome by considering only wing scatterings (Adams 1972). Hence, we are now writing a 
Monte Carlo code similar to that of Adams. Also we will investigate effects of following various 
deviations from our simplified configuration: the initial line width, the bulk motion of scatterers, 
the geometrical shape of scattering media, the degree of homogeneity of the dusts, and the 
distribution form of photon sources. 
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One of main goals of our research is to develop a tool to estimate the dust abundance in 
accordance with the usual method using the Balmer decrement (see Osterbrock 1989) and the 
slope of UV continuum (Fall & Pei 1989). This may help us correct properly the dust extinction 
of emission lines in starburst galaxies. Especially this will enable us to estimate the star formation 
rate of those galaxies which are located near or far. 

The other ramification is the formation of the damped Lya absorption (DLA) in many quasar 
spectra. The dust contents in DLA galaxies may deform the absorption profile. This may spoil 
the usual Voigt profile fitting procedure, and subsequent derivation of physical values can be 
erroneous. 

HML and SHA are grateful for the financial support of Brain Korea 21 project. HML was 
also supported by KOSEF grant no. 1999-2-113-001-5. HWL gratefully acknowledges support 
from the BK21 project. 
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